% COZZI AND IMPULLITTI (2015), GLOBALIZATION AND WAGE POLARIZATION, REVIEW
% OF ECONOMICS AND STATISTICS

% FUNCTION FILE FOR COMPUTATION OF EQUILIBRIUM AND COMPARATIVE STATICS


function f=funz(x,AF)


global sig1 gam beta phi AD tau epsi fi1 n lambda eta rho eserv fi zd1 zf1

   
the0D = x(1); % Ability cutoff for chosing education
the0F = x(2);
wDL = x(3);   % Unskilled wage D

wDH = (sig1*(the0D/(the0D-gam))*x(3)); % Skilled wage D
wFL = (((the0F-gam)/(sig1*the0F))^(1-beta)); % Unskilled wage F
wFH = (sig1*(the0F/(the0F-gam))*wFL);  %Skilled wage F

% ========================= TECHNOLOGY  ===================================

   
ZD = ((1/zd1)*(wDL^beta)*(wDH^(1-beta))); % Production cost function D
ZF = ((1/zf1)*(wFL^beta)*(wFH^(1-beta))); % Production cost function F
FD = ((1/AD)*((wDL^phi)*(wDH)^(1-phi)));  % Innovation cost function D
FF = ((1/AF)*((wFL^phi)*(wFH)^(1-phi)));  % Innovation cost function F

   % Marginal products: production
   
ZDL = ((1/zd1)*(beta*((wDH/wDL)^(1-beta))));
ZDH = ((1/zd1)*((1-beta)*((wDH/wDL)^(-beta))));
ZFL = ((1/zf1)*(beta*((wFH/wFL)^(1-beta))));
ZFH = (1/zf1)*((1-beta)*((wFH/wFL)^(-beta)));

   % Marginal products: R&D
   
FDL = ((phi*(1/AD)*((wDH/wDL)^(1-phi))));
FDH = (((1-phi)*(1/AD)*((wDH/wDL)^(-phi))));
FFL = ((phi*(1/AF)*((wFH/wFL)^(1-phi))));
FFH = (((1-phi)*(1/AF)*((wFH/wFL)^(-phi))));

% =========================== CUTOFFS =====================================

   % Skill cutoff to hire service workers
   
theHSD_fun = @(y) 1-y^epsi-(((wDH/wDL)*(1-eserv)*(y-gam))^epsi);
theHSF_fun = @(y) 1-y^epsi-(((wFH/wFL)*(1-eserv)*(y-gam))^epsi);

theHSD0 = 0.6;
theHSF0 = 0.6;

theHSD = csolve(theHSD_fun,theHSD0,[],1e-7,300);
theHSF = csolve(theHSF_fun,theHSF0,[],1e-7,300);

   % Skill cutoff for unskilled in production

theLSD = ((1-(theHSD^epsi))^(1/epsi));
theLSF = ((1-(theHSF^epsi))^(1/epsi));
  
   % Labor supplies
   
LDm = ((epsi/(1+epsi))*(the0D^(1+epsi)-theLSD^(1+epsi)));  % unskilled supply of production workers
LFm = ((epsi/(1+epsi))*(the0F^(1+epsi)-theLSF^(1+epsi)));

HD = ((1-the0D^epsi)*fi*(((epsi*((1-theHSD^(1+epsi))*(2 - eserv) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam)); % Skilled labor supply
HF = ((1-the0F^epsi)*fi*(((epsi*((1-theHSF^(1+epsi))*(2 - eserv) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam));

   % Average workers skills
   
theDm_avg = ((epsi/(1+epsi)) *((the0D)^(1+epsi)-(theLSD)^(1+epsi))/(the0D^epsi-theLSD^epsi)); %production
theFm_avg = ((epsi/(1+epsi)) *((the0F)^(1+epsi)-(theLSF)^(1+epsi))/(the0F^epsi-theLSF^epsi));


theDS_avg = ((epsi/(1+epsi)) *theLSD);  % Unskilled (services)
theFS_avg = ((epsi/(1+epsi)) *theLSF);

theDH_avg = (((epsi*((1-theHSD^(1+epsi)) + theHSD^(1+epsi) - the0D^(1+epsi)))/((1+epsi)*(1-the0D^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply
theFH_avg = (((epsi*((1-theHSF^(1+epsi)) + theHSF^(1+epsi) - the0F^(1+epsi)))/((1+epsi)*(1-the0F^epsi)))-gam);  % This average corrects for the increased supply of labour by skills: it is the average over the individuals of the averages per unit time, which means that those who work 2-eserv hours have their wages divided by (2-eserv). This is why we do not see (2-eserv) there, which of course we see it in the aggregate skill supply

theDunskill_avg = ((epsi/(1+epsi))*the0D); % Average unskilled (service+production)
theFunskill_avg = ((epsi/(1+epsi))*the0F);

% =================== OTHER ENDOGENOUS VARIABLES ==========================

% Here we report some expressions we use the reduced the original
% equilibrium system to 3 equations.

qF_til = ((1-fi1)*((lambda^(eta-1))-1)/n);
Fr = (FD/FF);
cD_til =( (Fr*(lambda-1)-(lambda-ZD*tau)) / ((lambda-ZD/tau)-(lambda-1)*Fr) );
x_til = ((1/FF)* (((lambda-1)/lambda)*(1+cD_til))/ (rho+(n/((1-fi1)*((lambda^(eta-1))-1))) + (n*fi1/(1-fi1))));
IF = (1/qF_til)*(1/(cD_til+1));
CF = ((1/IF)*(LFm / (qF_til*(ZFL/lambda)*(1+cD_til) + x_til*FFL)));
qF = (qF_til*IF);
qD = (1-(qF));
ID = ((1/qF_til)-IF);
CD = (qD/qF)*CF;
xindex = ((x_til)*CF); % R&D Difficulty indiex

% ===================== EQUILIBRIUM SYSTEM ================================


f(1) = HD- qD*(ZDH/lambda)*(CD/tau + CF*tau)-ID*xindex*FDH; % market clearing skilled D
f(2) = LDm- qD*(ZDL/lambda)*(CD/tau +CF*tau)-ID*xindex*FDL; % market clearing unkilled D
f(3) = HF - (LFm*(ZFH*qF_til*(1+cD_til)/lambda + x_til*FFH) / ( qF_til*((ZFL/lambda)*(1+cD_til))+x_til*FFL ));












